02_Normalize_Reads

Normalize reads and create DGE objects

library(tidyverse)
library(edgeR)

get reads and create data frame

files <- dir("../input/Kallisto_output/", include.dirs = TRUE)
files %>% head()

counts.list <- map(files, ~ read_tsv(
  file=file.path("..","input","Kallisto_output",.,"abundance.tsv"),
  col_types = "cdddd"))
names(counts.list) <- files

counts <- sapply(counts.list, select, est_counts) %>% 
  bind_cols(counts.list[[1]][,"target_id"],.)
counts[is.na(counts)] <- 0
colnames(counts) <- sub(".est_counts","",colnames(counts),fixed = TRUE)
counts

write_csv(counts,"../output/2018-timecourse_V3.0_raw_counts_.csv.gz")

create sample description data frame

key <- readxl::read_excel("../input/tube_no_legend_time_course_2018.xlsx",
                          na=c("","na"),
                          col_types=c("text", "text", "text", "skip", "text", "skip", "skip", "skip", "text", "text", "text", "skip", "text", "skip", "skip", "date", "date"))  %>%
  mutate(sampling_time_specific=format(sampling_time_specific, format="%H:%M:%S"))
key

create reformatted tube_no

key <- key %>%
  mutate(tube_no_2 = {
    tolower(tube_no) %>%
      str_replace("q_([1-9](_|$))", "q_00\\1") %>%
      str_replace("q_([1-9][0-9](_|$))", "q_0\\1") }) %>%
  select(tube_no, tube_no_2, everything())
key

samples <- tibble(
  file=files,
  tube_no_2 = str_extract(files, pattern = "q_[0-9]{3}(_d8)?") 
)
samples

samples <- left_join(samples, key) 
samples <- samples %>%
  mutate(sampling_day=str_pad(sampling_day,width=2,side = "left",pad="0")) %>%
  mutate(group=str_c(genotype, soil_trt, sampling_day, sampling_time,sep="-"))
samples

counts2 <- counts %>% 
  as.data.frame() %>% 
  column_to_rownames(var = "target_id") %>%
  as.matrix() %>%
  round(0)

samples2 <- samples %>%
  select(-tube_no_2, -tube_no, -pot, -plant_no, -sampling_day_specific, -sampling_time_specific) %>%
  as.data.frame() %>%
  column_to_rownames(var="file")

dge <- DGEList(counts=counts2, samples=samples2, group=samples2$group)

normalize

dge <- calcNormFactors(dge)

save(dge, file="../output/timecourseDGE.Rdata")

save some formatted cpm files

cpm averaged for each sample type

dim(log2cpmGroup)

[1] 46250    48

write.csv(log2cpmGroup, "../output/log2cpmGroup.csv.gz")

cpm for each individual sample

write.tsv(log2cpmSample, "../output/log2cpmSample.tsv.gz")

Error in write.tsv(log2cpmSample, "../output/log2cpmSample.tsv.gz") : 
  could not find function "write.tsv"